X=imread('grey.jpg');
x=double(X);
[C,L]=size(x);

%算子选择对话框
str = {'Sobel','Prewitt','Roberts','Marr','Canny'};
[type,ok] = listdlg('PromptString','Select a type:',...
                      'SelectionMode','single',...
                      'ListString',str,'ListSize', [200, 75]);
if ok == 0
    return ;
end

%手动计算出来的7×7的高度对称高斯滤波模板，主要用于canny算子的计算
div = 12.2791;
G = [    0.0111    0.0388    0.0821    0.1054    0.0821    0.0388    0.0111;
         0.0388    0.1353    0.2865    0.3679    0.2865    0.1353    0.0388;
         0.0821    0.2865    0.6065    0.7788    0.6065    0.2865    0.0821;
         0.1054    0.3679    0.7788    1.0000    0.7788    0.3679    0.1054;
         0.0821    0.2865    0.6065    0.7788    0.6065    0.2865    0.0821;
         0.0388    0.1353    0.2865    0.3679    0.2865    0.1353    0.0388;
         0.0111    0.0388    0.0821    0.1054    0.0821    0.0388    0.0111];

%必要的变量声明
result1 = 0;
result2 = 0;%一个记录横向梯度大小，一个记录纵向梯度大小
out = zeros(C,L);


switch type,
%%%%%%%%%%%%%%%%
  case 1,      %% Sobel算子
      factor =6;
      limit = 20;
      for I = 2:C-1
         for J = 2:L-1
            result1 = x(I+1,J-1)+2*x(I+1, J)+x(I+1,J+1)...%利用差分模板计算梯度
                      -x(I-1, J-1)-2*x(I-1,J)-x(I-1,J+1);
            result2 = x(I-1,J+1)+2*x(I, J+1)+x(I+1,J+1)...
                     -x(I-1, J-1)-2*x(I,J-1)-x(I+1,J-1);
            result = max(result1, result2);
            out(I,J) = result/factor;
            if out(I,J) <1
                out(I,J) = 1;
            end
            if out(I,J) >256
                out(I,J) = 256;
            end
         end
     end
     %二值化
     for I = 1:C
         for J = 1:L
              if out(I,J)>=limit
                  out(I,J) = 256;
              else
                  out(I,J) = 1;
              end
         end
     end
%%%%%%%%%%%%%%%%
  case 2,    %%   Prewitt算子
        factor = 3;
        limit = 25;
        for I = 2:C-1
           for J = 2:L-1 
               result1 = x(I+1,J-1)+x(I+1, J)+x(I+1,J+1)...%利用差分模板计算梯度
                   -x(I-1, J-1)-x(I-1,J)-x(I-1,J+1);
               result2 = x(I-1,J+1)+2*x(I, J+1)+x(I+1,J+1)...
                   -x(I-1, J-1)-2*x(I,J-1)-x(I+1,J-1); 
               result = max(result1, result2);
               out(I,J) = result/factor;
               if out(I,J) <1
                  out(I,J) = 1;
               end
               if out(I,J) >256
                  out(I,J) = 256;
               end
           end
        end
        %二值化
        for I = 1:C
           for J = 1:L
                if out(I,J)>=limit
                    out(I,J) = 256;
                else
                    out(I,J) = 1;
                end
           end
        end
%%%%%%%%%%%%%%%%
  case 3,     %%  Roberts算子
         factor = 1;
         limit =30 ;
         for I = 1:C-1
            for J = 2:L-1 
                result1 = x(I,J)-x(I+1,J-1); %利用差分模板计算梯度
                result2 = x(I,J)-x(I+1,J+1);
                result = max(result1, result2);
                out(I,J) = result/factor;
                if out(I,J) <1
                    out(I,J) = 1;
                end
                if out(I,J) >256
                  out(I,J) = 256;
                end
            end
         end
         %二值化
         for I = 1:C
           for J = 1:L
                if out(I,J)>=limit
                    out(I,J) = 256;
                else
                    out(I,J) = 1;
                end
           end
         end
%%%%%%%%%%%%%%%%
  case 4,            %%Marr
        factor = 16;
        limit = 4;
        temp = ones(C, L);
        %%%%%%%%%%%%%%%
        %%高斯滤波在本次实验中不是必要的，因为二阶求导模板已经包含高斯滤波成分。
        for I = 4:C-3
          for J = 4:L-3 
              result = 0;
              for u = -3:3
                 for v = -3:3
                     result = result + x(I+u, J+v)*G(u+4,v+4);
                 end
              end
              x(I,J) = result/div;
          end
        end
        %%%%%%%%%%%%%%%%
        %利用拉普拉斯高斯模板进行二阶求导
        for I = 3:C-2
          for J = 3:L-2 
              result = 16*x(I,J)-(x(I,J+1)+x(I,J-1)+x(I-1, J)+x(I+1, J))*2 ...
                    -x(I,J-2)-x(I,J+2)-x(I-1, J-1)-x(I-1, J+1)-x(I-2, J)...
                       -x(I+1,J-1)-x(I+1,J+1)-x(I-2, J);
               temp(I,J) = round(result/(factor*limit));
              end
        end
        %过零点检测，取边缘点
        for I = 3:C-2
           for J = 3:L-2
                    if temp(I, J-1)*temp(I, J+1) < 0 ||...
                            temp(I-1, J)*temp(I+1, J) < 0  ||...
                            temp(I-1, J-1)*temp(I+1, J+1) < 0 ||...
                            temp(I-1, J+1)*temp(I+1, J-1) < 0
                        out(I, J) = 256;
                    end
           end
        end
%%%%%%%%%%%%%%%%
   case 5,                %% Cany算子
        factor = 2;
        limit = 3.5;
        range = zeros(C, L);   %幅度矩阵
        tanangle = zeros(C, L); %角度矩阵(此处采用角度，也可直接用正切来比较)
        %%%%%%%%%%%%%%%
        %%高斯滤波
        for I = 4:C-3
          for J = 4:L-3 
              result = 0;
              for u = -3:3
                 for v = -3:3
                     result = result + x(I+u, J+v)*G(u+4,v+4);
                 end
              end
              x(I,J) = result/div;
          end
        end
        %%一阶差分卷积模板
        for I = 1:C-1
            for J = 1:L-1
                result1 = x(I+1,J)+x(I+1,J+1)-x(I,J)-x(I,J+1);
                result2 = x(I,J)+x(I+1,J)-x(I,J+1)-x(I+1,J+1);
                result = (result1*result1+result2*result2)^0.5;
                result=result/factor;
                range(I, J) = round(result);    %求出梯度大小   
                if result1 == 0
                    result1 = 0.000001;
                end
                tanangle(I, J) = atan(result2/result1);  %求出梯度的夹角
            end
        end
        
        %%非极大值抑制
        % 在[-pi/2, pi/2]上讨论，分四个方向
        rangetemp = range;
        for I = 2:C-1
            for J = 2:L-1
                if tanangle(I, J) >= -pi/8 && tanangle(I, J) < pi/8 && ( rangetemp(I, J)...
                        <=  rangetemp(I+1, J) || rangetemp(I, J) <=  rangetemp(I-1, J))
                    range(I, J) = 0;
                end
                if tanangle(I, J) >= pi/8 && tanangle(I, J) < 3*pi/8 && ( rangetemp(I, J)...
                        <=  rangetemp(I-1, J+1) ||  rangetemp(I, J) <=  rangetemp(I+1, J-1))
                    range(I, J) = 0;
                end
                if tanangle(I, J) >= -3*pi/8 && tanangle(I, J) < -pi/8 && ( rangetemp(I, J)...
                        <=  rangetemp(I-1, J-1) ||  rangetemp(I, J) <=  rangetemp(I+1, J+1))
                    range(I, J) = 0;
                end
                if (tanangle(I, J) < -3*pi/8 || tanangle(I, J) >=  3*pi/8) && ( rangetemp(I, J)...
                        <=  rangetemp(I-1, J) ||  rangetemp(I, J) <=  rangetemp(I+1, J))
                    range(I, J) = 0;
                end
            end
        end
        %%阈值化，第一次取边缘点
        for I = 2:C-1
            for J = 2:L-1
                if range(I, J) > limit*2
                    out(I, J) = 256;
                end
            end
        end
        %做100次扫描确保将边缘都连接起来，第二次取边缘点
        for K = 1:100
         for I = 2:C-1
            for J = 2:L-1
                if range(I, J) > limit && range(I, J) <= limit*2
                    if  out(I-1, J) == 256 || out(I-1, J-1) == 256 ||...    %边缘连接
                            out(I-1, J+1) == 256 ||  out(I+1, J-1) == 256 ||...
                             out(I+1, J) == 256 ||  out(I+1, J+1) == 256 ||...
                            out(I, J-1) == 256 ||  out(I, J+1) == 256 
                    
                        out(I, J) = 256;
                    end
                end
            end
         end
        end
%%%%%%%%%%%%%%%%
    otherwise,
end  



out = uint8(out)-1;
subplot(1,1,1);
imshow(out);
str = 'the edge of the image';
switch type,
    case 1,
       imwrite(out,'sobel边缘检测.jpg');
       str = 'the edge of the image with "sobel"'; 
    case 2,
       imwrite(out,'Prewitt边缘检测.jpg');
       str = 'the edge of the image with "Prewitt"'; 
    case 3,
       imwrite(out,'Roberts边缘检测.jpg');
       str = 'the edge of the image with "Roberts"'; 
    case 4,
       imwrite(out,'Marr边缘检测.jpg');
       str = 'the edge of the image with "Marr"';
    case 5,
       imwrite(out,'Canny边缘检测.jpg');
       str = 'the edge of the image with "Canny"'; 
end
title(str);

